#for flowchart
library(DiagrammeR)
#to export flowchart as .png
#webshot::install_phantomjs()
library(DiagrammeRsvg)
library(rsvg)
library(grid)
library(xtable)
library(dplyr)
library(ggplot2)
library(gridExtra)
#nest/unnest
library(tidyr)
#map function (kind of like a for loop)
library(purrr)
#tidy model summary
library(broom)
library(readxl)
#for tableby
library(arsenal)

#Also uses functions from plyr, scales, mgcv, cowplot

Define geom_split_violin()

Based on https://debruine.github.io/posts/plot-comparison/

GeomSplitViolin <- ggproto("GeomSplitViolin", GeomViolin, draw_group = function(self,
    data, ..., draw_quantiles = NULL) {
    data <- transform(data, xminv = x - violinwidth * (x - xmin), xmaxv = x + violinwidth *
        (xmax - x))
    grp <- data[1, "group"]
    newdata <- plyr::arrange(transform(data, x = if (grp%%2 == 1)
        xminv else xmaxv), if (grp%%2 == 1)
        y else -y)
    newdata <- rbind(newdata[1, ], newdata, newdata[nrow(newdata), ], newdata[1,
        ])
    newdata[c(1, nrow(newdata) - 1, nrow(newdata)), "x"] <- round(newdata[1, "x"])
    if (length(draw_quantiles) > 0 & !scales::zero_range(range(data$y))) {
        stopifnot(all(draw_quantiles >= 0), all(draw_quantiles <= 1))
        quantiles <- ggplot2:::create_quantile_segment_frame(data, draw_quantiles)
        aesthetics <- data[rep(1, nrow(quantiles)), setdiff(names(data), c("x", "y")),
            drop = FALSE]
        aesthetics$alpha <- rep(1, nrow(quantiles))
        both <- cbind(quantiles, aesthetics)
        quantile_grob <- GeomPath$draw_panel(both, ...)
        ggplot2:::ggname("geom_split_violin", grid::grobTree(GeomPolygon$draw_panel(newdata,
            ...), quantile_grob))
    } else {
        ggplot2:::ggname("geom_split_violin", GeomPolygon$draw_panel(newdata, ...))
    }
})

geom_split_violin <- function(mapping = NULL, data = NULL, stat = "ydensity", position = "identity",
    ..., draw_quantiles = NULL, trim = TRUE, scale = "area", na.rm = FALSE, show.legend = NA,
    inherit.aes = TRUE) {
    layer(data = data, mapping = mapping, stat = stat, geom = GeomSplitViolin, position = position,
        show.legend = show.legend, inherit.aes = inherit.aes, params = list(trim = trim,
            scale = scale, draw_quantiles = draw_quantiles, na.rm = na.rm, ...))
}

Load Data

load('../Data/DataWithPropensities_seed1.RData')

#convert PrimaryDiagnosis to factor
dat3$PrimaryDiagnosis <- factor(dat3$PrimaryDiagnosis, levels = c("Autism", "None"))
levels(dat3$PrimaryDiagnosis)
## [1] "Autism" "None"
dat3$PrimaryDiagnosis <- relevel(dat3$PrimaryDiagnosis, "None")
levels(dat3$PrimaryDiagnosis) = c("TD", "ASD")

49 ASD-only + 353 TD = 402 ADHD_Secondary=0

Reshape data to combine motion QC levels

# create dummy factor to include all subjects
dat3$noExclusion <- ifelse(dat3$ID > 0, "Pass", "Pass")
dat3$noExclusion <- factor(dat3$noExclusion, levels = c("Pass", "Fail"))

# convert KKI_criteria to factor with reference level 'Pass'
dat3$KKI_criteria <- factor(dat3$KKI_criteria, levels = c("Pass", "Fail"))

# convert Ciric_length to factor with reference level 'Pass' to match
# KKI_criteria
dat3$Ciric_length <- factor(dat3$Ciric_length, levels = c("Pass", "Fail"))

# combine Ciric_length, KKI, and noExclusion exclusion into one variable
allVariables = c("ID", "PrimaryDiagnosis", "AgeAtScan", "Ciric_length", "KKI_criteria",
    "noExclusion", "PANESS.TotalOverflowNotAccountingForAge", "SRS.Score", "WISC.GAI",
    "DuPaulHome.InattentionRaw", "DuPaulHome.HyperactivityRaw", "ADOS.Comparable.Total",
    "CurrentlyOnStimulants", "HeadCoil", "Sex", "ADHD_Secondary", "SES.Family", "Race2",
    "handedness", "CompletePredictorCases", "YearOfScan", "MeanFramewiseDisplacement.KKI")

idVariables = c("ID", "PrimaryDiagnosis", "AgeAtScan", "PANESS.TotalOverflowNotAccountingForAge",
    "SRS.Score", "WISC.GAI", "DuPaulHome.InattentionRaw", "DuPaulHome.HyperactivityRaw",
    "ADOS.Comparable.Total", "CurrentlyOnStimulants", "HeadCoil", "Sex", "ADHD_Secondary",
    "SES.Family", "Race2", "handedness", "CompletePredictorCases", "YearOfScan",
    "MeanFramewiseDisplacement.KKI")

qcMelt <- reshape2::melt(dat3[, allVariables], id.vars = names(dat3)[which(names(dat3) %in%
    idVariables)], variable.name = "Motion.Exclusion.Level", value.name = "Included")

# rename exclusion levels NOTE: need None to be highest level for
# geom_split_violin
levels(qcMelt$Motion.Exclusion.Level) <- c("Strict", "Lenient", "None")

# convert Included to factor with pass as reference
qcMelt$Included <- factor(qcMelt$Included, levels = c("Pass", "Fail"))

# rename levels of value
levels(qcMelt$Included) <- c("Included", "Excluded")

with(qcMelt, table(PrimaryDiagnosis, Motion.Exclusion.Level, Included))
## , , Included = Included
## 
##                 Motion.Exclusion.Level
## PrimaryDiagnosis Strict Lenient None
##              TD     151     308  372
##              ASD     29     114  173
## 
## , , Included = Excluded
## 
##                 Motion.Exclusion.Level
## PrimaryDiagnosis Strict Lenient None
##              TD     221      64    0
##              ASD    144      59    0

Lenient motion QC = KKI_criteria

Strict motion QC = Ciric_length

Limit initial data set to complete cases

dat3 <- filter(dat3, CompletePredictorCases==1)

1. Impact of motion QC on sample size

completeCases <- filter(qcMelt, CompletePredictorCases==1)

tabN<- tableby(Motion.Exclusion.Level~
                 Included, data=filter(completeCases, Motion.Exclusion.Level!="None"))
summary(tabN,  
        title='Proportion of complete cases included/excluded by motion QC',digits=0, digits.p=4,digits.pct=1, numeric.simplify=TRUE, total=FALSE, test=FALSE)
Proportion of complete cases included/excluded by motion QC
Strict (N=504) Lenient (N=504)
Included
   Included 171 (33.9%) 403 (80.0%)
   Excluded 333 (66.1%) 101 (20.0%)
tabASD<- tableby(Motion.Exclusion.Level~Included, data=filter(completeCases,
                                                                    PrimaryDiagnosis=="ASD"))

summary(tabASD,  
        title = "Proportion of ASD complete cases included/excluded",
        digits=0, digits.p=4,digits.pct=1, numeric.simplify=TRUE, total=FALSE, test=FALSE)
Proportion of ASD complete cases included/excluded
Strict (N=151) Lenient (N=151) None (N=151)
Included
   Included 29 (19.2%) 107 (70.9%) 151 (100.0%)
   Excluded 122 (80.8%) 44 (29.1%) 0 (0.0%)
tabTD<- tableby(Motion.Exclusion.Level~Included, data=filter(completeCases,
                                                                    PrimaryDiagnosis=="TD"))

summary(tabTD,  
        title = "Proportion of TD complete cases included/excluded",
        digits=0, digits.p=4,digits.pct=1, numeric.simplify=TRUE, total=FALSE, test=FALSE)
Proportion of TD complete cases included/excluded
Strict (N=353) Lenient (N=353) None (N=353)
Included
   Included 142 (40.2%) 296 (83.9%) 353 (100.0%)
   Excluded 211 (59.8%) 57 (16.1%) 0 (0.0%)

1b. Summarize sociodemographic info for complete cases for paper

#make M reference level for sex
completeCases$Sex <- relevel(as.factor(completeCases$Sex), "M")

#use chisq for Sex and handedness, kruskal-wallis rank test for Age
tabSex<- tableby( PrimaryDiagnosis~Sex+AgeAtScan+handedness, data=filter(completeCases, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(numeric.test="kwt", cat.test="chisq", total=FALSE))

#use fisher's exact test for Race, k-w for SES
tabRace<- tableby( PrimaryDiagnosis~Race2+SES.Family, data=filter(completeCases, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(numeric.test="kwt", cat.test="fe", total=FALSE))

tab12 <- merge(tabSex, tabRace)

#Currently on Stimulants - no test because no TDs are currently on stimulants by design
completeCases$CurrentlyOnStimulants <- factor(completeCases$CurrentlyOnStimulants)

#rename factor levels
levels(completeCases$CurrentlyOnStimulants)[levels(completeCases$CurrentlyOnStimulants)=="0"] <- "No"
levels(completeCases$CurrentlyOnStimulants)[levels(completeCases$CurrentlyOnStimulants)=="1"] <- "Yes"

tabStim<- tableby( PrimaryDiagnosis~CurrentlyOnStimulants, data=filter(completeCases, Motion.Exclusion.Level=="None" & Included=="Included"), control=tableby.control(total=FALSE, test=FALSE))

tab123 <- merge(tab12, tabStim)
summary(tab123)
TD (N=353) ASD (N=151) p value
Sex < 0.001
   M 245 (69.4%) 127 (84.1%)
   F 108 (30.6%) 24 (15.9%)
AgeAtScan 0.826
   Mean (SD) 10.363 (1.248) 10.324 (1.363)
   Range 8.020 - 12.980 8.010 - 12.990
handedness 0.253
   Left 17 (4.8%) 12 (7.9%)
   Mixed 19 (5.4%) 11 (7.3%)
   Right 317 (89.8%) 128 (84.8%)
Race2 0.004
   African American 36 (10.2%) 9 (6.0%)
   Asian 27 (7.6%) 3 (2.0%)
   Biracial 45 (12.7%) 12 (7.9%)
   Caucasian 245 (69.4%) 127 (84.1%)
SES.Family 0.006
   Mean (SD) 54.135 (9.390) 51.964 (9.379)
   Range 18.500 - 66.000 27.000 - 66.000
CurrentlyOnStimulants
   No 353 (100.0%) 97 (64.2%)
   Yes 0 (0.0%) 54 (35.8%)

Generate version of table to paste into overleaf

tab <- xtable(as.data.frame(summary(tab123)))
print(tab, type="latex")
## % latex table generated in R 4.1.1 by xtable 1.8-4 package
## % Thu Oct  7 00:43:14 2021
## \begin{table}[ht]
## \centering
## \begin{tabular}{rllll}
##   \hline
##  &  & TD (N=353) & ASD (N=151) & p value \\ 
##   \hline
## 1 & **Sex** &  &  & $<$ 0.001 \\ 
##   2 & \&nbsp;\&nbsp;\&nbsp;M & 245 (69.4\%) & 127 (84.1\%) &  \\ 
##   3 & \&nbsp;\&nbsp;\&nbsp;F & 108 (30.6\%) & 24 (15.9\%) &  \\ 
##   4 & **AgeAtScan** &  &  & 0.826 \\ 
##   5 & \&nbsp;\&nbsp;\&nbsp;Mean (SD) & 10.363 (1.248) & 10.324 (1.363) &  \\ 
##   6 & \&nbsp;\&nbsp;\&nbsp;Range & 8.020 - 12.980 & 8.010 - 12.990 &  \\ 
##   7 & **handedness** &  &  & 0.253 \\ 
##   8 & \&nbsp;\&nbsp;\&nbsp;Left & 17 (4.8\%) & 12 (7.9\%) &  \\ 
##   9 & \&nbsp;\&nbsp;\&nbsp;Mixed & 19 (5.4\%) & 11 (7.3\%) &  \\ 
##   10 & \&nbsp;\&nbsp;\&nbsp;Right & 317 (89.8\%) & 128 (84.8\%) &  \\ 
##   11 & **Race2** &  &  & 0.004 \\ 
##   12 & \&nbsp;\&nbsp;\&nbsp;African American & 36 (10.2\%) & 9 (6.0\%) &  \\ 
##   13 & \&nbsp;\&nbsp;\&nbsp;Asian & 27 (7.6\%) & 3 (2.0\%) &  \\ 
##   14 & \&nbsp;\&nbsp;\&nbsp;Biracial & 45 (12.7\%) & 12 (7.9\%) &  \\ 
##   15 & \&nbsp;\&nbsp;\&nbsp;Caucasian & 245 (69.4\%) & 127 (84.1\%) &  \\ 
##   16 & **SES.Family** &  &  & 0.006 \\ 
##   17 & \&nbsp;\&nbsp;\&nbsp;Mean (SD) & 54.135 (9.390) & 51.964 (9.379) &  \\ 
##   18 & \&nbsp;\&nbsp;\&nbsp;Range & 18.500 - 66.000 & 27.000 - 66.000 &  \\ 
##   19 & **CurrentlyOnStimulants** &  &  &  \\ 
##   20 & \&nbsp;\&nbsp;\&nbsp;No & 353 (100.0\%) & 97 (64.2\%) &  \\ 
##   21 & \&nbsp;\&nbsp;\&nbsp;Yes & 0 (0.0\%) & 54 (35.8\%) &  \\ 
##    \hline
## \end{tabular}
## \end{table}

1c. Create exclusion flowchart

1d. Proportion excluded

Define theme for proportion excluded plots

#Easier to make three separate panels and combine them than to use faceting function
My_Theme_prop = theme_light()+theme(
  legend.title =element_blank(),
  axis.title.x = element_text(size = 12),
  axis.title.y = element_text(size = 10),
  plot.title = element_text(size = 30),
  axis.text.x = element_text(size = 10),
  axis.text.y = element_text(size = 10),
  strip.text.x = element_text(size = 10,color="black"),
  strip.background = element_rect(fill="white"))

1d.1 Proportion excluded in each Diagnostic Group

motion <- filter(completeCases, Motion.Exclusion.Level != "None")
motion$Motion.Exclusion.Level <- droplevels(motion$Motion.Exclusion.Level)

motion <- group_by(motion, PrimaryDiagnosis, Motion.Exclusion.Level, Included)

dx_proportions <- ggplot(motion, aes(x = PrimaryDiagnosis, fill = Included)) + geom_bar(position = "fill",
    alpha = 0.6) + facet_grid(~Motion.Exclusion.Level) + scale_fill_manual(values = c("#FDE599",
    "#9FB0CC")) + scale_color_manual(values = c("#E9D38D", "#8C9AB4")) + My_Theme_prop +
    theme(legend.title = element_blank()) + theme(legend.title = element_blank()) +
    ylab("Proportion of Children") + theme(legend.position = "right") + theme(legend.margin = margin(t = 0,
    r = 0, b = -1, l = -1)) + theme(legend.key.size = unit(0.08, "in")) + theme(axis.title.x = element_blank()) +
    theme(axis.text.x = element_text(size = 10))

png("Figures/fig_propExcludedDx_cc.png", width = 3, height = 2, units = "in", res = 200)
dx_proportions
invisible(dev.off())

dx_proportions

dx_proportions <- dx_proportions + theme(legend.position = "none") + theme(axis.text.x = element_text(size = 8),
    axis.title.x = element_text(size = 12)) + xlab("Diagnosis Group")

# Pearson's chi squared tests
extib <- tibble(motion)

exNest <- extib %>%
    select(c("PrimaryDiagnosis", "Motion.Exclusion.Level", "Included")) %>%
    group_by(Motion.Exclusion.Level) %>%
    tidyr::nest()

# nested models
ex_chisq <- exNest %>%
    mutate(stats = map(data, ~broom::tidy(chisq.test(.x$PrimaryDiagnosis, .x$Included)))) %>%
    unnest(stats)

ex_chisq
## # A tibble: 2 × 6
## # Groups:   Motion.Exclusion.Level [2]
##   Motion.Exclusion.Level data               statistic  p.value parameter method 
##   <fct>                  <list>                 <dbl>    <dbl>     <int> <chr>  
## 1 Strict                 <tibble [504 × 2]>      19.9  8.07e-6         1 Pearso…
## 2 Lenient                <tibble [504 × 2]>      10.3  1.30e-3         1 Pearso…

The proportion of children excluded differs across diagnostic groups using both the lenient and strict motion QC.

1d.2 Currently on Stimulants

stim <- filter(motion, PrimaryDiagnosis=="ASD")
stim$PrimaryDiagnosis <- droplevels(stim$PrimaryDiagnosis)
stim$CurrentlyOnStimulants <- factor(stim$CurrentlyOnStimulants)
levels(stim$CurrentlyOnStimulants)
## [1] "No"  "Yes"
levels(stim$CurrentlyOnStimulants) = c("No", "Yes")

stimNest <- tibble(stim) %>% 
  select(c("CurrentlyOnStimulants",
           "Motion.Exclusion.Level", "Included")) %>% 
  group_by(Motion.Exclusion.Level) %>% 
  tidyr::nest()

#nested models
stim_chisq <- stimNest %>% 
  mutate(stats = map(data, ~broom::tidy(chisq.test(.x$CurrentlyOnStimulants, .x$Included)))) %>%
  unnest(stats)

stim_chisq
## # A tibble: 2 × 6
## # Groups:   Motion.Exclusion.Level [2]
##   Motion.Exclusion.Level data               statistic p.value parameter method  
##   <fct>                  <list>                 <dbl>   <dbl>     <int> <chr>   
## 1 Strict                 <tibble [151 × 2]>      6.40  0.0114         1 Pearson…
## 2 Lenient                <tibble [151 × 2]>      6.39  0.0115         1 Pearson…

1d.5 Proportion excluded by Currently on stimulants in ASD group

stim_proportions <- ggplot(stim, aes(x = CurrentlyOnStimulants, fill = Included))+
  geom_bar(position = "fill", alpha = .6)+
  facet_grid(~Motion.Exclusion.Level)+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme_prop+
  xlab("Currently On Stimulants?")+
  theme(legend.title = element_blank())+
  #theme(legend.position="top")+
  theme(legend.text = element_text(size=7))+
  ylab("Proportion of Autistic Children")+
  theme(legend.position = "none")+
  theme(strip.text.x = element_blank())
  #theme(plot.margin = unit(c(1, 6, 1, 6), "cm"))

prop_legend = cowplot::get_legend(stim_proportions + guides(color = guide_legend(ncol = 1))+
                                    theme(legend.position = "bottom", 
                                          legend.text = element_text(size = 11),
                                          legend.key.size=unit(.1, "in")))

1d.6 Combine proportion plot with legend & save

####Impact of motion QC on framewise displacement metrics

#mean and max FD are different for lenient and strict because some frames at the beginning and/or end of the scan are excluded for scans to pass lenient motion QC
dat3$MeanFD.None = dat3$MeanFramewiseDisplacement
dat3$MaxFD.None = dat3$MaxFramewiseDisplacement

#same for all levels
dat3$FramesWithFDLessThanOrEqualTo250microns.None = dat3$FramesWithFDLessThanOrEqualTo250microns
dat3$FramesWithFDLessThanOrEqualTo250microns.KKI = dat3$FramesWithFDLessThanOrEqualTo250microns

meanFD = c("ID", "MeanFramewiseDisplacement", "MeanFramewiseDisplacement.KKI", 
          "MeanFD.None")

maxFD = c("ID",  "MaxFramewiseDisplacement", "MaxFramewiseDisplacement.KKI", "MaxFD.None")

framesFD = c("ID",  "FramesWithFDLessThanOrEqualTo250microns",
             "FramesWithFDLessThanOrEqualTo250microns.KKI",
             "FramesWithFDLessThanOrEqualTo250microns.None")

fdID = c("ID")

meanFD.df <- reshape2::melt(dat3[, meanFD],
                         id.vars=names(dat3)[which(names(dat3) %in%  fdID)], 
                         variable.name = "Motion.Exclusion.Level",
                         value.name = "MeanFramewiseDisplacement")

levels(meanFD.df$Motion.Exclusion.Level)
## [1] "MeanFramewiseDisplacement"     "MeanFramewiseDisplacement.KKI"
## [3] "MeanFD.None"
#rename levels to match motion QC levels in completeCases
levels(meanFD.df$Motion.Exclusion.Level) <- c("Strict", "Lenient", "None")

#repeat for MaxFD
maxFD.df <- reshape2::melt(dat3[, maxFD],
                         id.vars=names(dat3)[which(names(dat3) %in%  fdID)], 
                         variable.name = "Motion.Exclusion.Level",
                         value.name = "MaxFramewiseDisplacement")

levels(maxFD.df$Motion.Exclusion.Level)
## [1] "MaxFramewiseDisplacement"     "MaxFramewiseDisplacement.KKI"
## [3] "MaxFD.None"
#rename levels to match motion QC levels in completeCases
levels(maxFD.df$Motion.Exclusion.Level) <- c("Strict", "Lenient", "None")

#merge meanFD.df and maxFD.df
fdMerg <- merge(meanFD.df, maxFD.df)

#repeat for FramesWithFDLessThanOrEqualTo250microns
frames.df <- reshape2::melt(dat3[, framesFD],
                         id.vars=names(dat3)[which(names(dat3) %in%  fdID)], 
                         variable.name = "Motion.Exclusion.Level",
                         value.name = "FramesWithFDLessThanOrEqualTo250microns")

levels(frames.df$Motion.Exclusion.Level)
## [1] "FramesWithFDLessThanOrEqualTo250microns"     
## [2] "FramesWithFDLessThanOrEqualTo250microns.KKI" 
## [3] "FramesWithFDLessThanOrEqualTo250microns.None"
#rename levels to match motion QC levels in qcMelt
levels(frames.df$Motion.Exclusion.Level) <- c("Strict", "Lenient", "None")

#merge with fdMerg
fdMerg <- merge(fdMerg, frames.df)

#merge with completeCases
completeCases <- merge(completeCases, fdMerg)

passOnly <- filter(completeCases, Included=="Included")

mean framewise displacement split violin none, lenient, strict

meanFD_violin <- ggplot(passOnly, aes(Motion.Exclusion.Level, MeanFramewiseDisplacement,
                                      fill = PrimaryDiagnosis, color=PrimaryDiagnosis)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 2) +
  stat_summary(fun = "mean", position = position_dodge(width = 0.18), 
               color="black", geom="point", aes(y=MeanFramewiseDisplacement))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.18), 
               color="black", geom="errorbar", width=.15)+
  scale_fill_manual(values = c("#009E73", "#FDE599"))+
  scale_color_manual(values = c("#05634a", "#E9D38D"))+
  My_Theme_prop+theme(legend.title =element_blank())+
  theme(axis.text.y=element_text(size=8))+
  theme(legend.text = element_text(size = 7))+
  theme(axis.title.y = element_blank())+
  theme(axis.title.x = element_blank())+
  theme(legend.position = c(0.35, .7))+
  theme(legend.key.size=unit(.15, "in"))+
  theme(legend.box.margin = margin(-2, -2, -2, -2))+
  ggtitle("Mean FD")+
  theme(plot.title = element_text(size = 8, hjust=0.5))

max framewise displacement split violin none, lenient, strict

maxFD_violin <- ggplot(passOnly, aes(Motion.Exclusion.Level, MaxFramewiseDisplacement,
                                     fill = PrimaryDiagnosis, color=PrimaryDiagnosis)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 2) +
  stat_summary(fun = "mean", position = position_dodge(width = 0.18), color="black", 
               geom="point", aes(y=MaxFramewiseDisplacement))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.18), 
               color="black", geom="errorbar", width=.15)+
  scale_fill_manual(values = c("#009E73", "#FDE599"))+
  scale_color_manual(values = c("#05634a", "#E9D38D"))+
  My_Theme_prop+
  theme(axis.text.y=element_text(size=8))+
  theme(axis.title.y = element_blank())+
  theme(axis.title.x = element_blank())+
  theme(legend.position = "none")+
  ggtitle("Max FD")+
  theme(plot.title = element_text(size = 8, hjust=0.5))

frames with FD <= 0.25 mm split violin none, lenient, strict

frames_violin <- ggplot(passOnly, aes(Motion.Exclusion.Level, FramesWithFDLessThanOrEqualTo250microns,
                                     fill = PrimaryDiagnosis, color=PrimaryDiagnosis)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 2) +
  stat_summary(fun = "mean", position = position_dodge(width = 0.18), color="black", 
               geom="point", aes(y=FramesWithFDLessThanOrEqualTo250microns))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.18), 
               color="black", geom="errorbar", width=.15)+
  scale_fill_manual(values = c("#009E73", "#FDE599"))+
  scale_color_manual(values = c("#05634a", "#E9D38D"))+
  My_Theme_prop+
  theme(axis.text.y=element_text(size=8))+
  theme(axis.title.y = element_blank())+
  theme(axis.title.x = element_blank())+
  theme(legend.position = "none")+
  ggtitle("Frames with FD<0.25 mm")+
  theme(plot.title = element_text(size = 8, hjust=0.5))

Combine 3 FD plots into one figure and save

## quartz_off_screen 
##                 2

Framewise displacement metrics are more similar across diagnostic groups using the strict motion QC, but very few participants are labeled as usable.

Filter out “None” from Motion.Exclusion.Level to make remaining group differences in FD metrics following motion QC easier to see

mean framewise displacement

passOnly <- filter(passOnly, Motion.Exclusion.Level!="None")
meanFD_violin <- ggplot(passOnly, aes(Motion.Exclusion.Level, MeanFramewiseDisplacement,
                                      fill = PrimaryDiagnosis, color=PrimaryDiagnosis)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 2) +
  stat_summary(fun = "mean", position = position_dodge(width = 0.18), 
               color="black", geom="point", aes(y=MeanFramewiseDisplacement))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.18), 
               color="black", geom="errorbar", width=.15)+
  scale_fill_manual(values = c("#009E73", "#FDE599"))+
  scale_color_manual(values = c("#05634a", "#E9D38D"))+
  My_Theme_prop+theme(legend.title =element_blank())+
  theme(axis.text.y=element_text(size=8))+
  theme(legend.text = element_text(size = 7))+
  theme(axis.title.y = element_blank())+
  theme(axis.title.x = element_blank())+
  theme(legend.position = c(0.35, .7))+
  theme(legend.key.size=unit(.15, "in"))+
  theme(legend.box.margin = margin(-2, -2, -2, -2))+
  ggtitle("Mean FD")+
  theme(plot.title = element_text(size = 8, hjust=0.5))

max framewise displacement

maxFD_violin <- ggplot(passOnly, aes(Motion.Exclusion.Level, MaxFramewiseDisplacement,
                                     fill = PrimaryDiagnosis, color=PrimaryDiagnosis)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 2) +
  stat_summary(fun = "mean", position = position_dodge(width = 0.18), color="black", 
               geom="point", aes(y=MaxFramewiseDisplacement))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.18), 
               color="black", geom="errorbar", width=.15)+
  scale_fill_manual(values = c("#009E73", "#FDE599"))+
  scale_color_manual(values = c("#05634a", "#E9D38D"))+
  My_Theme_prop+
  theme(axis.text.y=element_text(size=8))+
  theme(axis.title.y = element_blank())+
  theme(axis.title.x = element_blank())+
  theme(legend.position = "none")+
  ggtitle("Max FD")+
  theme(plot.title = element_text(size = 8, hjust=0.5))

frames with FD <= 0.25 mm

frames_violin <- ggplot(passOnly, aes(Motion.Exclusion.Level, FramesWithFDLessThanOrEqualTo250microns,
                                     fill = PrimaryDiagnosis, color=PrimaryDiagnosis)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 2) +
  stat_summary(fun = "mean", position = position_dodge(width = 0.18), color="black", 
               geom="point", aes(y=FramesWithFDLessThanOrEqualTo250microns))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.18), 
               color="black", geom="errorbar", width=.15)+
  scale_fill_manual(values = c("#009E73", "#FDE599"))+
  scale_color_manual(values = c("#05634a", "#E9D38D"))+
  My_Theme_prop+
  theme(axis.text.y=element_text(size=8))+
  theme(axis.title.y = element_blank())+
  theme(axis.title.x = element_blank())+
  theme(legend.position = "none")+
  ggtitle("Frames with FD<0.25 mm")+
  theme(plot.title = element_text(size = 8, hjust=0.5))

Combine 3 FD plots with just lenient and strict Motion.Exclusion.Levels and save

## quartz_off_screen 
##                 2

2. rs-fMRI exclusion probability as a function of age and symptom severity

Covariates

phenoVariables <- c("ID", "PrimaryDiagnosis", 
                    "ADOS.Comparable.Total", 
                    "SRS.Score", 
                    "PANESS.TotalOverflowNotAccountingForAge", 
                    "DuPaulHome.InattentionRaw",
                    "DuPaulHome.HyperactivityRaw", 
                    "AgeAtScan", 
                    "WISC.GAI",
                    "Motion.Exclusion.Level", "Included")

phenoIDs <- c("ID", "PrimaryDiagnosis", "Motion.Exclusion.Level", "Included")

aim1 <- reshape2::melt(completeCases[, phenoVariables],
                         id.vars=names(completeCases)[which(names(completeCases) %in% phenoIDs)])

levels(aim1$variable) <- c("ADOS", "SRS", "Motor Overflow", "Inattention", 
                           "Hyperactivity", "Age", "GAI")

aim1G <- group_by(aim1, PrimaryDiagnosis, Motion.Exclusion.Level, Included, variable)

2a. Univarate GAMs

aim1$delta = rep(NA,length=nrow(aim1))
aim1$delta = ifelse(aim1$Included=="Included",1,0)
aim1tib <- tibble(filter(aim1, Motion.Exclusion.Level!="None"))
aim1tib$Motion.Exclusion.Level <- droplevels(aim1tib$Motion.Exclusion.Level)

aim1Nest <- aim1tib %>% 
  group_by(Motion.Exclusion.Level, variable) %>% 
  tidyr::nest()

#nested models
nested_gams <- aim1Nest %>% 
  mutate(model = map(data, ~mgcv::gam(1-delta~s(value, k=-1), data = na.omit(.x), 
                                      family=binomial(link=logit), method="REML")), 
         coefs = map(model, tidy, conf.int = FALSE),
         Rsq = map_dbl(model, ~summary(.)$r.sq)) %>% 
  unnest(coefs)

#Ben: correct for 7 lenient and 7 strict 
nested_gams_len <-  nested_gams %>% 
  filter(Motion.Exclusion.Level=="Lenient")

nested_gams_len$p.fdr = p.adjust(nested_gams_len$p.value, method = "BH")

nested_gams_strict <-  nested_gams %>% 
  filter(Motion.Exclusion.Level=="Strict")
  
nested_gams_strict$p.fdr = p.adjust(nested_gams_strict$p.value, method = "BH")

#combine to print
nested_gams <- rbind(nested_gams_len, nested_gams_strict)

#list adjusted p values
nested_gams[, c(1:2,6:11)]
## # A tibble: 14 × 8
## # Groups:   Motion.Exclusion.Level, variable [14]
##    Motion.Exclusion.… variable     edf ref.df statistic  p.value    Rsq    p.fdr
##    <fct>              <fct>      <dbl>  <dbl>     <dbl>    <dbl>  <dbl>    <dbl>
##  1 Lenient            ADOS        2.04   2.47     18.4   2.08e-4 0.0442  4.86e-4
##  2 Lenient            SRS         1.83   2.29     17.4   2.95e-4 0.0466  5.17e-4
##  3 Lenient            Motor Ove…  1.00   1.00     15.6   7.58e-5 0.0308  4.54e-4
##  4 Lenient            Inattenti…  1.38   1.67      7.44  1.17e-2 0.0157  1.17e-2
##  5 Lenient            Hyperacti…  2.15   2.70     13.2   4.59e-3 0.0242  5.36e-3
##  6 Lenient            Age         1.00   1.00      9.33  2.26e-3 0.0164  3.17e-3
##  7 Lenient            GAI         1.00   1.00     14.7   1.30e-4 0.0288  4.54e-4
##  8 Strict             ADOS        1.00   1.00     21.9   2.79e-6 0.0444  9.76e-6
##  9 Strict             SRS         1.00   1.00     22.7   1.53e-6 0.0567  9.76e-6
## 10 Strict             Motor Ove…  1.00   1.00     10.2   1.42e-3 0.0194  1.99e-3
## 11 Strict             Inattenti…  1.00   1.00     17.8   2.47e-5 0.0360  4.31e-5
## 12 Strict             Hyperacti…  1.88   2.35     25.0   1.30e-5 0.0516  3.04e-5
## 13 Strict             Age         1.76   2.20      7.73  2.84e-2 0.0129  2.84e-2
## 14 Strict             GAI         1.00   1.00      7.55  6.02e-3 0.0130  7.02e-3
#max p value for 7 lenient models
max(nested_gams_len$p.fdr)
## [1] 0.01171046
#max p value for 7 strict models
max(nested_gams_strict$p.fdr)
## [1] 0.02839014
nested_gams <- nested_gams %>% 
           mutate(LB = map(data, ~round(min(na.omit(.x$value)))),
                  UB = map(data, ~round(max(na.omit(.x$value)))),
                  range = map2(LB, UB, ~seq(from=.x, to=.y, by=1)),
                  logpredict = map2(model, range, ~predict(.x, newdata = data.frame(value = .y), type="link",se.fit=TRUE)),
                  fit = map(logpredict, ~plogis(.x$fit)),
                  lCI = map(logpredict, ~plogis(.x$fit-1.96*.x$se.fit)),
                  hCI = map(logpredict, ~plogis(.x$fit+1.96*.x$se.fit)))

Define theme for gam plots

2a.1 Probability of exclusion as a function of ADOS

ados <- nested_gams %>% 
  filter(variable=="ADOS") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_ados <- ggplot(ados, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='Probability of Exclusion', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("ADOS")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))

2a.2 Probability of exclusion as a function of SRS

srs <- nested_gams %>% 
  filter(variable=="SRS") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_srs <- ggplot(srs, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+
  gam_theme+
  ggtitle("SRS")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.3 Probability of exclusion as a function of Inattention

ina <- nested_gams %>% 
  filter(variable=="Inattention") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_in <- ggplot(ina, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("Inattention")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.4 Probability of exclusion as a function of Hyperactivity

hi <- nested_gams %>% 
  filter(variable=="Hyperactivity") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_hi <- ggplot(hi, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("Hyperactivity")+
theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.5 Probability of exclusion as a function of Motor Overflow

mo <- nested_gams %>% 
  filter(variable=="Motor Overflow") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_mo <- ggplot(mo, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("Motor Overflow")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.6 Probability of exclusion as a function of Age

age<- nested_gams %>% 
  filter(variable=="Age") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_age <- ggplot(age, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("Age")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

2a.7 Probability of exclusion as a function of GAI

gai <- nested_gams %>% 
  filter(variable=="GAI") %>% 
  select("variable", "Motion.Exclusion.Level", "range", "fit", 'lCI', 'hCI') %>% 
  unnest(c(range, fit, lCI, hCI))

p_gai <- ggplot(gai, aes(x=range, y=fit))+
  geom_line(aes(colour = Motion.Exclusion.Level),size=1.2)+ylim(0,1)+theme_bw()+
  geom_ribbon(aes(ymin=lCI, ymax=hCI, fill=Motion.Exclusion.Level), linetype='blank', alpha=0.2)+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154","#9FB0CC"))+
  labs(x='', y='', fill='Motion Control', colour='Motion Control')+  
  scale_x_continuous(expand = c(0, 0))+ 
  gam_theme+
  ggtitle("GAI")+
  theme(plot.title = element_text(size = 11, hjust = 0.5))+
  theme(axis.title.y = element_blank())

p_legend = cowplot::get_legend(p_gai + guides(color = guide_legend(nrow = 1))+theme(legend.position = "bottom", legend.text = element_text(size = 10),   legend.key.size=unit(.1, "in")))

2b. Density plots of covariates used to fit GAMs (across included & excluded children)

Define theme for density plots of covariates across included and excluded children

2b.1 ADOS density

ddata <- nested_gams %>% 
  filter(variable=="ADOS") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data) %>% 
  filter(PrimaryDiagnosis=="ASD")

ddata$PrimaryDiagnosis <- droplevels(ddata$PrimaryDiagnosis)

d_ados=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0), limits = c(0, .09), breaks=seq(0, .08, by=.02))+  
  labs(x='', y='Density')+
  scale_fill_manual(values = c("#FDE599"))+ 
  scale_color_manual(values = c("#E9D38D"))+ 
  den_theme

2b.2 SRS density

ddata <- nested_gams %>% 
  filter(variable=="SRS") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_srs=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0), limits=c(0,max(srs$range)),breaks = seq(0, 100 , by = 50))+  
  scale_y_continuous(expand = c(0, 0))+  
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.3 Inattention density

ddata <- nested_gams %>% 
  filter(variable=="Inattention") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_in=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  labs(x='')+  
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  theme(axis.title.y = element_blank())+
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.4 Hyperactivity/Impulsivity Density

ddata <- nested_gams %>% 
  filter(variable=="Hyperactivity") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_hi=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.5 Motor Overflow Density

ddata <- nested_gams %>% 
  filter(variable=="Motor Overflow") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_mo=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  theme(axis.title.y = element_blank())+
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.6 Age Density

ddata <- nested_gams %>% 
  filter(variable=="Age") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_age=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0), limits=c(8,13), breaks = seq(8, 13 , by = 1))+  
  scale_y_continuous(expand = c(0, 0), limits=c(0,.29), breaks=seq(0, .25, by = .05))+ 
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  theme(axis.title.y = element_blank())+
  labs(x='')

2b.7 GAI Density

ddata <- nested_gams %>% 
  filter(variable=="GAI") %>% 
  filter(Motion.Exclusion.Level=="Lenient") %>% 
  select("variable", "Motion.Exclusion.Level", "data") %>% 
  unnest(data)

d_gai=ggplot(ddata, aes(x=value, fill=PrimaryDiagnosis, color=PrimaryDiagnosis))+  
  geom_density(alpha=0.5, inherit.aes=TRUE)+  
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0), limits = c(0, .035), breaks=seq(0., .03, by=.01))+  
  scale_fill_manual(labels=c('TD','ASD'), values = c("#009E73", "#FDE599"))+
  scale_color_manual(labels=c('TD','ASD'), values = c("#05634a", "#E9D38D"))+
  
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  den_theme+
  labs(x='')+  
  theme(axis.title.y = element_blank())

hist_legend = cowplot::get_legend(d_gai + guides(color = guide_legend(nrow = 1))+theme(legend.position = "bottom", legend.text = element_text(size = 10),   legend.key.size=unit(.1, "in")))

combine gam plots with densities & print

top_row <- cowplot::plot_grid(p_ados, p_srs, p_in, p_hi, p_mo, p_age, p_gai, ncol=7, rel_widths=c(1.18/7, .97/7, .97/7, .97/7, .97/7, .97/7, .97/7))
bottom_row <- cowplot::plot_grid(d_ados, d_srs, d_in, d_hi, d_mo, d_age, d_gai, ncol=7, rel_widths=c(1.18/7, .97/7, .97/7, .97/7, .97/7, .97/7, .97/7))
## Warning: Removed 91 rows containing non-finite values (stat_density).
## Warning: Removed 19 rows containing non-finite values (stat_density).
png("Figures/fig_probExclusion_allGAM_TD_ASD_cc.png",width=10,height=6,units="in",res=200)
cowplot::plot_grid(p_legend, top_row, NULL, bottom_row, NULL, hist_legend, nrow=6, rel_heights=c(.1, 1, -.07, .5, -.07, .1))
dev.off()
## quartz_off_screen 
##                 2
#png("~/Dropbox/Apps/Overleaf/MotionSelectionBias_rsfMRI/Figures/fig_probExclusion_allGAM_TD_ASD_cc.png",width=10,height=6,units="in",res=200)
cowplot::plot_grid(p_legend, top_row, NULL, bottom_row, NULL, hist_legend, nrow=6, rel_heights=c(.1, 1, -.07, .5, -.07, .1))

#dev.off()

NOTE: 19 rows = # of participants with missing Motor Overflow scores (imputed for deconfounded group difference)

3. Impact of motion QC on phenotypic representation within diagnostic groups

3a. Split violin plots

My_Theme = theme_light()+theme(
  legend.title = element_blank(),
  axis.title.x = element_blank(),
  axis.title.y = element_blank(),
  axis.text.x = element_text(size = 5),
  axis.text.y = element_text(size = 8),
  strip.text.x = element_text(size = 12, face = "bold", color="black"),
  strip.text.y = element_text(size = 10, color="black"),
  strip.background = element_rect(fill="white"),
  plot.title = element_text(size = 9, hjust = 0.5))

3a.1 ADOS (ASD only) split violin

NOTE: Missing ADOS scores for one participant evaluated at CARD

ados <- aim1 %>% 
  filter(PrimaryDiagnosis=="ASD" & variable=="ADOS") %>% 
  dplyr::select(-c("PrimaryDiagnosis", "variable"))

ados_violin <- ggplot(ados, aes(Motion.Exclusion.Level, value, fill = Included, color = Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE, adjust = 1.5) +
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(legend.text = element_text(size = 8))+
  theme(legend.position = "none")+
  ggtitle("ADOS")

3a.2 SRS split violin

srs <- filter(aim1G, variable=="SRS")
aim1_srs <- ggplot(srs, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~., scales = "fixed")+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  # geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y =element_blank())+
  theme(legend.position = "none")+
  ggtitle("SRS")
  
v_legend = cowplot::get_legend(aim1_srs + guides(color = guide_legend(nrow = 2))+
                                 theme(legend.position = "left", 
                                       legend.text = element_text(size = 10), 
                                       legend.key.size=unit(.1, "in")))
## Warning: Removed 273 rows containing non-finite values (stat_ydensity).
## Warning: Removed 273 rows containing non-finite values (stat_summary).

## Warning: Removed 273 rows containing non-finite values (stat_summary).

NOTE: 273/3 = 91, # of participants missing SRS

with(dat3, table(PrimaryDiagnosis, is.na(SRS.Score)))
##                 
## PrimaryDiagnosis FALSE TRUE
##              TD    263   90
##              ASD   150    1

3a.3 Inattention Split Violin

inat <- filter(aim1G, variable=="Inattention")
paper_inat <- ggplot(inat, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y = element_blank())+
  theme(legend.position = "none")+
  ggtitle("Inattention")

3a.4 Hyperactivity/Impulsivity Spilt Violin

hyp <- filter(aim1G, variable=="Hyperactivity")
paper_hyp <- ggplot(hyp, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y = element_blank())+
  theme(legend.position = "none")+
  ggtitle("Hyperactivity")

3a.5 Motor Overflow Split Violin

overflow <- filter(aim1G, variable=="Motor Overflow")
aim1_of <- ggplot(overflow, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y =element_blank())+
  theme(legend.position = "none")+
  ggtitle("Motor Overflow")

3a.6 Age Split Violin

age <- filter(aim1G, variable=="Age")
aim1_age <- ggplot(age, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(strip.text.y =element_blank())+
  theme(legend.position = "none")+
  ggtitle("Age")

3a.7 GAI Split Violin

gai <- filter(aim1G, variable=="GAI")
aim1_gai <- ggplot(gai, aes(Motion.Exclusion.Level, value, fill = Included, color=Included)) +
  geom_split_violin(trim=TRUE, alpha = 0.5, inherit.aes = TRUE) +
  facet_grid(PrimaryDiagnosis~.)+
  stat_summary(fun = "mean", position = position_dodge(width = 0.5), 
               color="black", geom="point", aes(y=value))+
  stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), 
               color="black", geom="errorbar", width=.2)+
  # geom_jitter(position=position_jitterdodge(jitter.height = .25), alpha = .3)+
  #geom_pointrange(
  #  data = summary_pheno,
  #  aes(Motion.Exclusion.Level, mean, ymin=min, ymax=max),
  #  shape = 20,
  #  position = position_dodge(width = 0.9))+
  scale_fill_manual(values = c("#FDE599", "#9FB0CC"))+
  scale_color_manual(values = c("#E9D38D", "#8C9AB4"))+
  My_Theme+
  theme(legend.position = "none")+
  ggtitle("GAI")

Combine split violins into one figure & save

3b. Mann-Whitney U tests to compare included vs excluded participants using lenient motion QC (13 tests)

#run lenient tests first
aim1tib <- tibble(filter(aim1, Motion.Exclusion.Level=="Lenient"))

aim1MW <- aim1tib %>% 
  group_by(PrimaryDiagnosis, variable) %>% 
  tidyr::nest()

#hypothesis: included children will have less severe symptoms. NOTE: ADOS only collected in ASD (9 tests)
nested_mw_less <- aim1MW %>% 
  filter(variable %in% c("SRS", "Inattention", "Hyperactivity",
  "Motor Overflow")|(variable=="ADOS"&PrimaryDiagnosis=="ASD")) %>% 
  mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="less", data = na.omit(.x))),
         includedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Included"))),
         excludedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Excluded"))),
  coefs = map(mwm, tidy)) %>% 
  unnest(c(coefs, includedMean, excludedMean))

#hypothesis: included children will be older and have greater GAI (4 tests)
nested_mw_greater<- aim1MW %>% 
  filter(variable %in% c("Age", "GAI")) %>% 
  mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="greater", data = na.omit(.x))),
         includedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Included"))),
         excludedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Excluded"))),
  coefs = map(mwm, tidy)) %>% 
  unnest(c(coefs, includedMean, excludedMean))

complete_mw <- rbind(nested_mw_less, nested_mw_greater)

complete_mw$p.fdr <- p.adjust(complete_mw$p.value, method = "BH")

complete_mw[order(complete_mw$PrimaryDiagnosis, decreasing=TRUE), c(1:2, 5:9)]
## # A tibble: 13 × 7
## # Groups:   PrimaryDiagnosis, variable [13]
##    PrimaryDiagnosis variable  includedMean excludedMean statistic p.value method
##    <fct>            <fct>            <dbl>        <dbl>     <dbl>   <dbl> <chr> 
##  1 ASD              ADOS             14.5         14.5      1780. 9.21e-3 Wilco…
##  2 ASD              SRS              NA           NA        1720  5.79e-3 Wilco…
##  3 ASD              Motor Ov…        NA           NA        1260  9.48e-4 Wilco…
##  4 ASD              Inattent…        16.7         16.7      2442. 6.42e-1 Wilco…
##  5 ASD              Hyperact…        12.0         12.0      2419  6.06e-1 Wilco…
##  6 ASD              Age              10.3         10.3      2848  2.16e-2 Wilco…
##  7 ASD              GAI             107.         107.       2960  6.55e-3 Wilco…
##  8 TD               SRS              NA           NA        4562. 4.31e-1 Wilco…
##  9 TD               Motor Ov…        NA           NA        7086. 5.69e-2 Wilco…
## 10 TD               Inattent…         3.10         3.10     8004  2.68e-1 Wilco…
## 11 TD               Hyperact…         2.13         2.13     7018. 1.98e-2 Wilco…
## 12 TD               Age              10.4         10.4     10074. 1.02e-2 Wilco…
## 13 TD               GAI             116.         116.       9882  2.02e-2 Wilco…

3c. Mann-Whitney U tests to compare included vs excluded participants using strict motion QC (13 tests)

#run tests using strict motion QC
aim1tib <- tibble(filter(aim1, Motion.Exclusion.Level=="Strict"))

aim1MW <- aim1tib %>% 
  group_by(variable, PrimaryDiagnosis) %>% 
  tidyr::nest()

#hypothesis: included children will have less severe symptoms
nested_mw_less <- aim1MW %>% 
  filter(variable %in% c("SRS", "Inattention", "Hyperactivity",
  "Motor Overflow")|(variable=="ADOS"&PrimaryDiagnosis=="ASD")) %>% 
  mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="less", data = na.omit(.x))),
         includedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Included"))),
         excludedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Excluded"))),
         coefs = map(mwm, tidy)) %>% 
  unnest(c(coefs, includedMean, excludedMean))

#hypothesis: included children will be older and have greater GAI
nested_mw_greater<- aim1MW %>% 
  filter(variable %in% c("Age", "GAI")) %>% 
  mutate(mwm = map(data, ~wilcox.test(value~Included, alternative="greater", data = na.omit(.x))),
         includedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Included"))),
         excludedMean = map(data, ~mean(.x$value, data=filter(.x, Included=="Excluded"))),
  coefs = map(mwm, tidy)) %>% 
  unnest(c(coefs, includedMean, excludedMean))

complete_mw <- rbind(nested_mw_less, nested_mw_greater)

complete_mw$p.fdr <- p.adjust(complete_mw$p.value, method = "BH")

complete_mw[order(complete_mw$PrimaryDiagnosis, decreasing=TRUE), c(1:2, 5:6, 10:11)]
## # A tibble: 13 × 6
## # Groups:   PrimaryDiagnosis, variable [13]
##    PrimaryDiagnosis variable       includedMean excludedMean alternative  p.fdr
##    <fct>            <fct>                 <dbl>        <dbl> <chr>        <dbl>
##  1 ASD              ADOS                  14.5         14.5  less        0.0767
##  2 ASD              SRS                   NA           NA    less        0.0763
##  3 ASD              Motor Overflow        NA           NA    less        0.0847
##  4 ASD              Inattention           16.7         16.7  less        0.223 
##  5 ASD              Hyperactivity         12.0         12.0  less        0.322 
##  6 ASD              Age                   10.3         10.3  greater     0.0763
##  7 ASD              GAI                  107.         107.   greater     0.303 
##  8 TD               SRS                   NA           NA    less        0.0847
##  9 TD               Motor Overflow        NA           NA    less        0.194 
## 10 TD               Inattention            3.10         3.10 less        0.194 
## 11 TD               Hyperactivity          2.13         2.13 less        0.0159
## 12 TD               Age                   10.4         10.4  greater     0.0847
## 13 TD               GAI                  116.         116.   greater     0.111
xtable(complete_mw[order(complete_mw$PrimaryDiagnosis, decreasing=TRUE), c(1:2, 5:6, 10:11)])
## % latex table generated in R 4.1.1 by xtable 1.8-4 package
## % Thu Oct  7 00:44:19 2021
## \begin{table}[ht]
## \centering
## \begin{tabular}{rllrrlr}
##   \hline
##  & PrimaryDiagnosis & variable & includedMean & excludedMean & alternative & p.fdr \\ 
##   \hline
## 1 & ASD & ADOS & 14.54 & 14.54 & less & 0.08 \\ 
##   2 & ASD & SRS &  &  & less & 0.08 \\ 
##   3 & ASD & Motor Overflow &  &  & less & 0.08 \\ 
##   4 & ASD & Inattention & 16.67 & 16.67 & less & 0.22 \\ 
##   5 & ASD & Hyperactivity & 11.99 & 11.99 & less & 0.32 \\ 
##   6 & ASD & Age & 10.32 & 10.32 & greater & 0.08 \\ 
##   7 & ASD & GAI & 107.11 & 107.11 & greater & 0.30 \\ 
##   8 & TD & SRS &  &  & less & 0.08 \\ 
##   9 & TD & Motor Overflow &  &  & less & 0.19 \\ 
##   10 & TD & Inattention & 3.10 & 3.10 & less & 0.19 \\ 
##   11 & TD & Hyperactivity & 2.13 & 2.13 & less & 0.02 \\ 
##   12 & TD & Age & 10.36 & 10.36 & greater & 0.08 \\ 
##   13 & TD & GAI & 115.54 & 115.54 & greater & 0.11 \\ 
##    \hline
## \end{tabular}
## \end{table}

4. Edgewise functional connectivity as a function of the covariates

Combine partial correlations with melted rs-fMRI usability and covariate info

startEdgeidx = which(names(dat3)=='r.ic1.ic2')
endEdgeidx = which(names(dat3)=='r.ic29.ic30')
names(dat3)[startEdgeidx:endEdgeidx]
##   [1] "r.ic1.ic2"   "r.ic1.ic4"   "r.ic1.ic8"   "r.ic1.ic13"  "r.ic1.ic14" 
##   [6] "r.ic1.ic15"  "r.ic1.ic17"  "r.ic1.ic19"  "r.ic1.ic21"  "r.ic1.ic22" 
##  [11] "r.ic1.ic24"  "r.ic1.ic25"  "r.ic1.ic26"  "r.ic1.ic27"  "r.ic1.ic28" 
##  [16] "r.ic1.ic29"  "r.ic1.ic30"  "r.ic2.ic4"   "r.ic2.ic8"   "r.ic2.ic13" 
##  [21] "r.ic2.ic14"  "r.ic2.ic15"  "r.ic2.ic17"  "r.ic2.ic19"  "r.ic2.ic21" 
##  [26] "r.ic2.ic22"  "r.ic2.ic24"  "r.ic2.ic25"  "r.ic2.ic26"  "r.ic2.ic27" 
##  [31] "r.ic2.ic28"  "r.ic2.ic29"  "r.ic2.ic30"  "r.ic4.ic8"   "r.ic4.ic13" 
##  [36] "r.ic4.ic14"  "r.ic4.ic15"  "r.ic4.ic17"  "r.ic4.ic19"  "r.ic4.ic21" 
##  [41] "r.ic4.ic22"  "r.ic4.ic24"  "r.ic4.ic25"  "r.ic4.ic26"  "r.ic4.ic27" 
##  [46] "r.ic4.ic28"  "r.ic4.ic29"  "r.ic4.ic30"  "r.ic8.ic13"  "r.ic8.ic14" 
##  [51] "r.ic8.ic15"  "r.ic8.ic17"  "r.ic8.ic19"  "r.ic8.ic21"  "r.ic8.ic22" 
##  [56] "r.ic8.ic24"  "r.ic8.ic25"  "r.ic8.ic26"  "r.ic8.ic27"  "r.ic8.ic28" 
##  [61] "r.ic8.ic29"  "r.ic8.ic30"  "r.ic13.ic14" "r.ic13.ic15" "r.ic13.ic17"
##  [66] "r.ic13.ic19" "r.ic13.ic21" "r.ic13.ic22" "r.ic13.ic24" "r.ic13.ic25"
##  [71] "r.ic13.ic26" "r.ic13.ic27" "r.ic13.ic28" "r.ic13.ic29" "r.ic13.ic30"
##  [76] "r.ic14.ic15" "r.ic14.ic17" "r.ic14.ic19" "r.ic14.ic21" "r.ic14.ic22"
##  [81] "r.ic14.ic24" "r.ic14.ic25" "r.ic14.ic26" "r.ic14.ic27" "r.ic14.ic28"
##  [86] "r.ic14.ic29" "r.ic14.ic30" "r.ic15.ic17" "r.ic15.ic19" "r.ic15.ic21"
##  [91] "r.ic15.ic22" "r.ic15.ic24" "r.ic15.ic25" "r.ic15.ic26" "r.ic15.ic27"
##  [96] "r.ic15.ic28" "r.ic15.ic29" "r.ic15.ic30" "r.ic17.ic19" "r.ic17.ic21"
## [101] "r.ic17.ic22" "r.ic17.ic24" "r.ic17.ic25" "r.ic17.ic26" "r.ic17.ic27"
## [106] "r.ic17.ic28" "r.ic17.ic29" "r.ic17.ic30" "r.ic19.ic21" "r.ic19.ic22"
## [111] "r.ic19.ic24" "r.ic19.ic25" "r.ic19.ic26" "r.ic19.ic27" "r.ic19.ic28"
## [116] "r.ic19.ic29" "r.ic19.ic30" "r.ic21.ic22" "r.ic21.ic24" "r.ic21.ic25"
## [121] "r.ic21.ic26" "r.ic21.ic27" "r.ic21.ic28" "r.ic21.ic29" "r.ic21.ic30"
## [126] "r.ic22.ic24" "r.ic22.ic25" "r.ic22.ic26" "r.ic22.ic27" "r.ic22.ic28"
## [131] "r.ic22.ic29" "r.ic22.ic30" "r.ic24.ic25" "r.ic24.ic26" "r.ic24.ic27"
## [136] "r.ic24.ic28" "r.ic24.ic29" "r.ic24.ic30" "r.ic25.ic26" "r.ic25.ic27"
## [141] "r.ic25.ic28" "r.ic25.ic29" "r.ic25.ic30" "r.ic26.ic27" "r.ic26.ic28"
## [146] "r.ic26.ic29" "r.ic26.ic30" "r.ic27.ic28" "r.ic27.ic29" "r.ic27.ic30"
## [151] "r.ic28.ic29" "r.ic28.ic30" "r.ic29.ic30"
signalFC <- dat3[, c(1,startEdgeidx:endEdgeidx)]

dat2 <- merge(aim1, signalFC, all=TRUE, by = "ID")
fcMelt <- reshape2::melt(dat2[,1:160],
                         id.vars=names(dat2)[1:7],
                         variable.name = "edge",
                         value.name = "fc")

4a. Run nested gams

fctib <- tibble(filter(fcMelt, Motion.Exclusion.Level!="None"))
fctib$Motion.Exclusion.Level <- droplevels(fctib$Motion.Exclusion.Level)

fcNest <- fctib %>% 
  filter(Included=="Included") %>% 
  group_by(variable, Motion.Exclusion.Level, edge) %>% 
  tidyr::nest()

#nested models
fc_gams <- fcNest %>% 
  mutate(model = map(data, ~mgcv::gam(fc~s(value), data = na.omit(.x), method="REML")), 
         coefs = map(model, tidy, conf.int = FALSE)) %>% 
  unnest(coefs)

4b. Plot histograms of edgewise p-values from GAMs

4b.1 Histogram of edgewise p-values for partial correlations as a function of ADOS.

NOTE: TD scores = 0

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="ADOS") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_ados_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  ggtitle("ADOS")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  labs(x='p value', y='Count')+
  theme(axis.title.y = element_text(size = 9, angle=90))+
  theme(axis.title.x = element_text(size = 7))

4b.2 Histogram of edgewise p-values for partial correlations as a function of SRS

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="SRS") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_srs_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("SRS")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.3 Histogram of p-values for partial correlations as a function of inattention

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="Inattention") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_in_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("Inattention")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.4 Histogram of p-values for partial correlations as a function of hyperactivity

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="Hyperactivity") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_hi_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("Hyperactivity")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.5 Histogram of p-values for partial correlations as a function of motor overflow

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="Motor Overflow") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_mo_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("Motor Overflow")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.6 Histogram of p-values for partial correlations as a function of age

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="Age") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_age_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("Age")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

4b.7 Histogram of p-values for partial correlations as a function of GAI

pdat <- fc_gams %>% 
  ungroup() %>% 
  filter(variable=="GAI") %>% 
  select(Motion.Exclusion.Level, p.value) %>% 
  group_by(Motion.Exclusion.Level)

hist_gai_p=ggplot(pdat, aes(x=p.value, fill=Motion.Exclusion.Level, color=Motion.Exclusion.Level))+  
  geom_histogram(position = "identity", alpha=0.5, inherit.aes=TRUE, binwidth = .05)+
  scale_x_continuous(expand = c(0, 0))+  
  scale_y_continuous(expand = c(0, 0))+  
  labs(x='p value', y='')+
  scale_color_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  scale_fill_manual(labels=c('Strict', 'Lenient'), values = c("#f55154", "#9FB0CC"))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(), panel.background = element_blank(), 
        axis.line = element_line(size=.5), panel.border = element_blank())+  
  My_Theme+
  theme(axis.title.y = element_blank())+
  ggtitle("GAI")+
  theme(plot.title = element_text(size = 9, hjust = 0.5))+
  theme(legend.position = "none")+
  theme(axis.title.x = element_text(size = 7))

hist_p_legend <- cowplot::get_legend(hist_gai_p + guides(color = guide_legend(nrow = 1))+theme(legend.position = "bottom", legend.text = element_text(size = 10),   legend.key.size=unit(.1, "in"), legend.title = element_blank()))

4b.8 Combine histograms and print

fc_hist <- cowplot::plot_grid(hist_ados_p, hist_srs_p, hist_in_p, hist_hi_p, hist_mo_p, hist_age_p, hist_gai_p, ncol=7, rel_widths=c(1.18/7, .97/7, .97/7, .97/7, .97/7, .97/7, .97/7))

png("Figures/fig_hist_rfc_cc.png",width=8,height=3,units="in",res=200)

cowplot::plot_grid(fc_hist, hist_p_legend, nrow=2, rel_heights=c(1, .1))
dev.off()
## quartz_off_screen 
##                 2
cowplot::plot_grid(fc_hist, hist_p_legend, nrow=2, rel_heights=c(1, .1))